function  [ rational_moments,  behavioral_moments ] = simulation_moments( sim, model_sol, plot_figure )
% Usage: moments_from_simulation( sim )
%
% This function generates key moments from simulated results from the model simulations.   The object sim is an output
% of the function "simulate.m".

par = model_sol.par;
dt = par.dt;

% Note: all growths are expressed as a yearly variable.

if(nargin<=2)
    plot_figure = false;
end

%%  ----------------------------------------------
%                                   Iterations
% ----------------------------------------------------
for(iter = 1:2)

if(iter==1)
    current_sim = sim.rational;
else
    if( model_sol.par.model_name=="behavioral" )  % if without the option to not have the behavioral case.
        current_sim = sim.behavioral;
    else
        behavioral_moments = {}  ;
        continue;
    end
end

%% Deal with NAN. This occurs sometimes because the growth in the economy is too fast (AH and AL too large) that the GDP at certain point is infinity.
output_growth_test = [current_sim.GDP(2:end)./current_sim.GDP(1:(end-1))-1;  zeros(1,1) ] ;
index_not_NA = find(~isnan(output_growth_test)& ~isinf(current_sim.GDP)) ;
if( mod(length(index_not_NA),1/dt)~=0  )
    index_not_NA =   index_not_NA( (mod(length(index_not_NA),1/dt)+1):end ) ;
end

current_sim.bank_equity = current_sim.bank_equity(index_not_NA);
current_sim.credit_spread = current_sim.credit_spread(index_not_NA);
current_sim.GDP = current_sim.GDP(index_not_NA);
current_sim.K = current_sim.K(index_not_NA);
current_sim.lambda = current_sim.lambda(index_not_NA);
current_sim.p = current_sim.p(index_not_NA);
current_sim.productivity = current_sim.productivity(index_not_NA);
current_sim.psi = current_sim.psi(index_not_NA);
current_sim.vol = current_sim.vol(index_not_NA);
current_sim.vol_premia = current_sim.vol_premia(index_not_NA);
current_sim.w = current_sim.w(index_not_NA);
current_sim.wb = current_sim.wb(index_not_NA);
current_sim.xK = current_sim.xK(index_not_NA);
current_sim.yK = current_sim.yK(index_not_NA);

N_sim = length(current_sim.GDP);
dN_effective = sim.dN_vec(index_not_NA);
N_interval = floor(1/par.dt);  % how many data points per year
dN_temp = movmean( dN_effective , [0, 11] );
dN_yearly = dN_temp(1:12:end)>0 ;


%% Basic Moments
output_growth_monthly_vec = [current_sim.GDP(2:end)./current_sim.GDP(1:(end-1))-1;  zeros(1,1) ] ;
output_growth_yearly_at_monthly_frequency =  [current_sim.GDP(13:end)./current_sim.GDP(1:(end-12))-1;  zeros(12,1) ] ;
GDP_sim = current_sim.GDP;
GDP_sim_yearly = GDP_sim(1:N_interval:end);
output_growth_yearly_vec = [GDP_sim_yearly(2:end)./GDP_sim_yearly(1:(end-1))-1; 0] ;
output_growth_3year_vec = [GDP_sim_yearly(4:end)./GDP_sim_yearly(1:(end-3))-1; zeros(3,1) ] ;
output_growth_avg = mean(output_growth_yearly_vec)   ;
output_growth_vol =  std(output_growth_yearly_vec);

% % Credit spread 
credit_spread_sim = current_sim.credit_spread ;  

% others
bank_equity_sim =  current_sim.w .* current_sim.p .* current_sim.K ;
bank_equity_return_sim =   [ bank_equity_sim(2:end)./bank_equity_sim(1:(end-1)) - 1  ;   0 ]*12 + par.rho ;    % Note: this is a annualized bank equity return, taking into the transition attrition and the dividend payment (consumption) by equity holders
equity_yearly = bank_equity_sim(1:N_interval:end) ;
bank_equity_excess_yearly_return_sim = [bank_equity_sim(13:end)./bank_equity_sim(1:(end-12)) - 1;  zeros(12,1) ] + movsum( (par.rho -  current_sim.equity_transition_returns - current_sim.rf)*dt , [0,12]) ;
bank_equity_excess_yearly_return_sim = bank_equity_excess_yearly_return_sim - mean(bank_equity_excess_yearly_return_sim) + mean(current_sim.equity_risk_premium_vec)*12;  % Adjustment because the directly calculated levels are much more accurate, depite the inaccuracy at fluctuations.
bank_equity_return_sim_yearly =  [equity_yearly(2:end)./equity_yearly(1:(end-1)) - 1;  0] + par.rho  -  current_sim.equity_transition_returns(1:N_interval:end);
bank_equity_excess_return_sim_yearly = [equity_yearly(2:end)./equity_yearly(1:(end-1)) - 1;  0] + par.rho -  current_sim.equity_transition_returns(1:N_interval:end) - current_sim.rf(1:12:end);
bank_equity_excess_return_sim_yearly = bank_equity_excess_return_sim_yearly - mean(bank_equity_excess_return_sim_yearly) + mean(current_sim.equity_risk_premium_vec)*12;  % Adjustment because the directly calculated levels are much more accurate, depite the inaccuracy at fluctuations.

bank_total_credit_sim = current_sim.w .* current_sim.p .* current_sim.K .* (current_sim.xK-1) ;    % bank credit = bank debt from our definition. 
bank_credit_ratio_sim =  bank_total_credit_sim  ./ current_sim.GDP ;   % bank credit to GDP ratio. Important!
bank_credit_ratio_sim = bank_credit_ratio_sim/std(bank_credit_ratio_sim);   % Normalize. 
bank_credit_ratio_growth_monthly = [ bank_credit_ratio_sim(2:end) -  bank_credit_ratio_sim(1:end-1) ;  0 ] ;
bank_credit_growth_monthly = [ bank_total_credit_sim(2:end) -  bank_total_credit_sim(1:end-1) ;  0 ] ./ bank_total_credit_sim ;

bank_credit_growth_sim_yearly = [bank_total_credit_sim(13:end)./bank_total_credit_sim(1:(end-12)) - 1;  0];
bank_credit_growth_sim_yearly = min(bank_credit_growth_sim_yearly(1:N_interval:end), quantile(bank_credit_growth_sim_yearly, 0.995));

bank_equity_yearly_return_avg = mean(bank_equity_return_sim_yearly);
bank_equity_yearly_return_vol = std(bank_equity_return_sim_yearly);
bank_leverage_sim = min( current_sim.xK, 20 );

%%  ----------------------------------------------
%                         Define a Crisis and Recessions
% ----------------------------------------------------
% Find the quantile that result in a yearly frequency of crisis at 4%. 
crisis_quantile = 0.04;  
x = bank_credit_ratio_growth_monthly;
% x = bank_credit_growth_monthly;
% x = output_growth_monthly_vec;

% *** Previous definition of crisis
% frequency_yearly_fun = @(q)  1000*  (sum( sum( reshape( (x<=quantile(x, q) ) , 1/par.dt , []) ) >0 ) / (N_sim*dt)  - crisis_quantile); 
% q_monthly = fzero( frequency_yearly_fun,  [0, crisis_quantile] );  % Note: using the interval is much better!!!!
% is_crisis_vec = x<=quantile(x, q_monthly);

% *** Another definition by removing the sparse crises. 
distance = 3*12;  % 3 years minimum distance between cirses. 
residual_fun = @(x_cutoff)  mean(gen_sparse_crisis(  x<x_cutoff ,  distance ))*12 - crisis_quantile;
x_cutoff = fzero( residual_fun,  quantile(x, crisis_quantile/10)  );
is_crisis_vec = gen_sparse_crisis(  x<x_cutoff ,  distance );
if( model_sol.par.model_name~="behavioral" | (model_sol.par.model_name=="behavioral"&&iter==2) )
    disp( ['sparse crises frequency = ',  num2str(100*crisis_quantile),  '%' ]  );
    disp( ['non-sparse crises frequency = ',  num2str(100*mean(x<x_cutoff)*12),  '%' ]  );
end
frequency_non_sparse_crisis = mean(x<x_cutoff)*12;

% **** define equity crashes 
equity_monthly_growth = [bank_equity_sim(2:end)./bank_equity_sim(1:(end-1))-1; 0]; % bank equity growth 
residual_fun = @(equity_growth_cutoff)  mean(gen_sparse_crisis(  equity_monthly_growth<equity_growth_cutoff ,  distance ))*12 - crisis_quantile;
equity_growth_cutoff = fzero( residual_fun,  quantile(equity_monthly_growth, crisis_quantile/10)  );
is_equity_crash_vec = gen_sparse_crisis(  equity_monthly_growth<equity_growth_cutoff ,  distance );
if( model_sol.par.model_name~="behavioral" | (model_sol.par.model_name=="behavioral"&&iter==2) )
    disp( ['sparse equity crash frequency = ',  num2str(100*crisis_quantile),  '%' ]  );
    disp( ['non-sparse equity crash frequency = ',  num2str(100*mean(equity_monthly_growth<equity_growth_cutoff)*12),  '%' ]  );
end
frequency_non_sparse_equity_crashes = mean(equity_monthly_growth<equity_growth_cutoff)*12;
% *** end of defining equity crahses

crises_starting_index_vec =  find(is_crisis_vec) ;
yearly_shocks = reshape( is_crisis_vec, 1/par.dt ,  []) ; % arbitrary data

% % % ------ Temporary Changes, explore the different bank credit ratio definitions  ------
% % quantile(bank_credit_ratio_growth_monthly, 0.0034)
% % quantile(output_growth_monthly_vec, 0.0034)
% plot( output_growth_monthly_vec(is_crisis_vec) , bank_credit_ratio_growth_monthly(is_crisis_vec), 'ro' );  xlabel('output growth'); ylabel('bank credit/GDP change'); hold on;
% plot( output_growth_monthly_vec(~is_crisis_vec) , bank_credit_ratio_growth_monthly(~is_crisis_vec), 'b+' );   
% legend( { 'Crisis',  'Non Crisis'  }  )
% % % ------- Temporary Changes ------

is_crisis_yearly_vec  =  (  sum(yearly_shocks)>0 )' ;
crises_starting_index_yearly_vec = find(is_crisis_yearly_vec);
% is_crisis_yearly_vec = output_growth_yearly_vec<=quantile(output_growth_yearly_vec, crisis_quantile);
% crises_starting_index_yearly_vec = find(is_crisis_yearly_vec);
yearly_distress_shocks = reshape( dN_effective, 1/par.dt ,  []) ; % arbitrary data
is_distress_yearly_vec = (  sum(yearly_distress_shocks)>0 )' ;

%  Average Frequency of Distress
frequency_crises_simulated = sum(is_crisis_yearly_vec) / length(is_crisis_yearly_vec) ;
frequency_distress_simulated = sum(dN_effective) / sim.T_vec(end) ;
ratio = par.lambdaHL / ( par.lambdaHL + par.lambdaLH );
frequency_distress_underlying = par.lambda_H * (1-ratio) + par.lambda_L * ratio;


% Recessions
% output_growth_next_year =  movsum( output_growth_monthly_vec ,  [0, 11] );
is_recession_yearly_vec = output_growth_yearly_vec < mean(output_growth_yearly_vec) - 0.01;
is_non_fin_recession_yearly_vec = is_recession_yearly_vec & ~is_distress_yearly_vec ; 
% Then convert them to numeric value
is_recession_yearly_vec = double(is_recession_yearly_vec);
is_non_fin_recession_yearly_vec = double(is_non_fin_recession_yearly_vec);
% is_recession_vec = zeros(size(output_growth_next_year)); 
% i=1;
% while( i < numel(is_recession_vec) )
%     if( output_growth_next_year(i) < mean(output_growth_next_year) - 0.015 )
%         is_recession_vec(i) = 1;  i = i + 12;
%     else
%         i = i+1;
%     end
% end
% sum(is_recession_vec) / length(is_recession_vec)  


%% TEST --- Comment before formal running
% index = 10000:11000; figure(1);
% output_growth_vec = [GDP_sim(2:end)./GDP_sim(1:(end-1))-1; 0] ;   % output drop in future 2 months
% productivity_growth_vec = [productivity_sim(2:end)./productivity_sim(1:(end-1))-1; 0];
% plot( current_sim.w(1000:3000), current_sim.lambda(1000:3000) )
% 
% figure; plot( current_sim.spread_basis(current_sim.spread_basis<0.9) , output_growth_vec(current_sim.spread_basis<0.9), 'o' )
% 
% figure; plot( current_sim.spread_basis(current_sim.spread_basis<0.05) , output_growth_vec(current_sim.spread_basis<0.05), 'o' )
% 
% % find the abnormal 
% index_abnormal = find( current_sim.spread_basis<0.05 &  output_growth_vec>-0.14 & output_growth_vec<-0.04 ); 
% w_abnormal = current_sim.w(index_abnormal); 
% lambda_abnormal = current_sim.lambda(index_abnormal); 
% dN_effective(index_abnormal)
% lambda_pre = current_sim.lambda(index_abnormal - par.lag_years/par.dt );
% 
% output_drop1 = model_sol.kappa_output_fitted_matrix( w_abnormal, lambda_abnormal )
% output_drop2 = model_sol.kappa_output_fitted_tensor( w_abnormal, lambda_abnormal, lambda_pre )
% output_drop_actual = output_growth_vec(index_abnormal)
% productivity_drop_actual = productivity_growth_vec(index_abnormal)
% 
% subindex = model_sol.spread_basis_matrix<0.05 &  model_sol.spread_basis_matrix>0.025 & model_sol.w_matrix<max(current_sim.w) ;
% figure; plot( model_sol.spread_basis_matrix(subindex) ,  -model_sol.kappa_output_matrix(subindex),   'o'  ); 
% model_sol.kappa_output_matrix(subindex)
% 
% 
% figure(2);
% subplot(1,2,1); surf(  w_matrix,  lambda_matrix, model_sol.kappa_output_matrix); xlabel('w'); ylabel('lambda'); zlabel('kappa output');  xlim( [0, 0.25] );
% subplot(1,2,2); surf(  w_matrix,  lambda_matrix, model_sol.spread_basis_matrix); xlabel('w'); ylabel('lambda'); zlabel('credit spread basic definition');  xlim( [0, 0.25] );
% 
% plot(w_grid, model_sol.kappa_output_matrix)
% plot(lambda_grid, model_sol.kappa_output_matrix(find(w_grid<max(current_sim.w)),:) )
% plot(lambda_grid, model_sol.spread_basis_matrix(find(w_grid<max(current_sim.w)),:) )
%% END OF TEST --- Comment before formal running


%% Output Declines During Crises in a 1 year horizon. Pick the middle point in the next year to proxy for the one-year average. 
output_change_3years_in_crises_vec=0; 
output_change_in_crises_vec=0;
bank_equity_change_in_crises_vec=0;
credit_spread_change_in_crises_vec=0;
for( year_index = 1:12 )
    output_change_3years_in_crises_vec = ( GDP_sim( min(crises_starting_index_vec+12*2+year_index, N_sim) ) ./ GDP_sim(crises_starting_index_vec) - 1 )/12 + output_change_3years_in_crises_vec; 
    output_change_in_crises_vec = ( GDP_sim( min(crises_starting_index_vec+year_index, N_sim) ) ./ GDP_sim(crises_starting_index_vec) - 1 )/12  + output_change_in_crises_vec ; 
    bank_equity_change_in_crises_vec = ( bank_equity_sim( min(crises_starting_index_vec+year_index, N_sim) ) ./ bank_equity_sim(crises_starting_index_vec)  - 1 )/12 +  bank_equity_change_in_crises_vec; 
    credit_spread_change_in_crises_vec = ( credit_spread_sim( min((crises_starting_index_vec+year_index), N_sim) ) - credit_spread_sim(crises_starting_index_vec) )/12 +  credit_spread_change_in_crises_vec; 
end
output_change_in_crises_avg = mean(output_change_in_crises_vec);   % In one year horizon. 
bank_equity_change_in_crises_avg = mean(bank_equity_change_in_crises_vec); 
credit_spread_change_in_crises_avg = mean(credit_spread_change_in_crises_vec); 

crises_trough_index_vec = crises_starting_index_vec;
full_recovery_index_vec = crises_starting_index_vec;
credit_spread_recovery_index_vec = crises_starting_index_vec;
GDP_recovery_index_vec = crises_starting_index_vec;
bank_equity_recovery_index_vec = crises_starting_index_vec;

% plot( T_vec, GDP_sim ,  T_vec, GDP_sim_smoothed)
avg_path_init = zeros(20*12, 1);
average_path.GDP = avg_path_init;
average_path.credit_spread = avg_path_init;
average_path.bank_equity = avg_path_init;
average_path.update_counts = avg_path_init;
average_path.lambda = avg_path_init;

prob_cutoff = 0.75;
for( i = 1:numel(crises_starting_index_vec) )
    start_index = crises_starting_index_vec(i);
%     plot( credit_spread_sim( start_index: (start_index+100) ) )

    next_start = crises_starting_index_vec(min(i+1,end));
    trough_index = start_index;
    recovery_index = start_index+1;
    credit_spread_recovery_index = start_index + 1; 
    GDP_recovery_index = start_index + 1; 
    bank_equity_recovery_index = start_index + 1; 
    while( trough_index<N_sim && GDP_sim(trough_index)> GDP_sim(trough_index+1) )  
        trough_index = trough_index + 1; 
    end
%     [~,lowest_index] = min(GDP_sim(start_index:min((start_index+8*N_interval),N_sim)));
%     trough_index = lowest_index-1 + start_index;

    while( recovery_index<N_sim && GDP_sim(recovery_index)< GDP_sim(start_index) && dN_effective(recovery_index)~=1  )  
        recovery_index = recovery_index + 1; 
    end
    credit_spread_change = credit_spread_sim(start_index+1) - credit_spread_sim(start_index);
    while(  credit_spread_recovery_index<N_sim && credit_spread_recovery_index< next_start && credit_spread_sim(credit_spread_recovery_index)> credit_spread_sim(start_index) +  credit_spread_change*prob_cutoff ) 
        credit_spread_recovery_index = credit_spread_recovery_index + 1; 
        if(credit_spread_recovery_index>= next_start ||  dN_effective(credit_spread_recovery_index)==1 )
            credit_spread_recovery_index = nan;
        end
    end
%     if( ~isnan(credit_spread_recovery_index)  )
%         figure;
%         years_vec = [1:( credit_spread_recovery_index + 84 - start_index+1 ) ]  * dt;
%         plot(  years_vec,  credit_spread_sim( start_index:(credit_spread_recovery_index+84) )    );
%         pause(2);
%     end
    
    GDP_change = GDP_sim(start_index+1) - GDP_sim(start_index);
    while(  GDP_recovery_index<N_sim && GDP_sim(GDP_recovery_index) < GDP_sim(start_index) +  GDP_change*prob_cutoff ) 
        GDP_recovery_index = GDP_recovery_index + 1; 
        if(GDP_recovery_index>= next_start || dN_effective(GDP_recovery_index)==1  )
            GDP_recovery_index = nan;
        end
    end
    equity_change = bank_equity_sim(start_index+1) - bank_equity_sim(start_index);
    while(  bank_equity_recovery_index<N_sim && bank_equity_recovery_index< next_start && bank_equity_sim(bank_equity_recovery_index) < bank_equity_sim(start_index) +  equity_change*prob_cutoff ) 
        bank_equity_recovery_index = bank_equity_recovery_index + 1; 
        if(bank_equity_recovery_index>= next_start  ||  dN_effective(bank_equity_recovery_index)==1  )
            bank_equity_recovery_index = nan;
        end
    end
    crises_trough_index_vec(i) = trough_index;
    full_recovery_index_vec(i) = recovery_index;
    credit_spread_recovery_index_vec(i) = credit_spread_recovery_index;
    GDP_recovery_index_vec(i) = GDP_recovery_index;
    bank_equity_recovery_index_vec(i) = bank_equity_recovery_index;
    
    % average path
    path_end_index = start_index + 1;
    while(  path_end_index< min(N_sim,next_start) && dN_effective(path_end_index)~=1  ) 
        path_end_index = path_end_index + 1; 
    end
    
    path_start_index = max(start_index - 24,1);  % a year ago before a crisis
    end_avg_path_index = min(  path_start_index+length(average_path.GDP)-1 , next_start );
    N_new = (end_avg_path_index-path_start_index+1);
    
    % Note: the average paths are in reference to the crisis starting year
    gdp_change_wrt_reference = current_sim.GDP(path_start_index:end_avg_path_index)./current_sim.GDP(start_index)-1;
    credit_spread_change_wrt_reference = current_sim.credit_spread(path_start_index:end_avg_path_index)-current_sim.credit_spread(start_index);
    bank_equity_change_wrt_reference = current_sim.bank_equity(path_start_index:end_avg_path_index)./current_sim.bank_equity(start_index)-1;
    lambda_chnage_wrt_reference =  current_sim.lambda(path_start_index:end_avg_path_index) - current_sim.lambda(start_index);
    
    average_path.GDP(1:N_new) = ( average_path.update_counts(1:N_new).*average_path.GDP(1:N_new) + gdp_change_wrt_reference ) ./ ( average_path.update_counts(1:N_new) + 1 ) ;
    average_path.credit_spread(1:N_new) = ( average_path.update_counts(1:N_new).*average_path.credit_spread(1:N_new) + credit_spread_change_wrt_reference ) ./ ( average_path.update_counts(1:N_new) + 1 ) ;
    average_path.bank_equity(1:N_new) = ( average_path.update_counts(1:N_new).*average_path.bank_equity(1:N_new) + bank_equity_change_wrt_reference ) ./ ( average_path.update_counts(1:N_new) + 1 ) ;
    average_path.lambda(1:N_new) =     ( average_path.update_counts(1:N_new).*average_path.lambda(1:N_new) + lambda_chnage_wrt_reference ) ./ ( average_path.update_counts(1:N_new) + 1 ) ;
    
    average_path.update_counts(1:N_new) = average_path.update_counts(1:N_new) + 1;
end
peak_to_trough_vec = (crises_trough_index_vec - crises_starting_index_vec) * dt; 
recovery_length_vec = (full_recovery_index_vec - crises_starting_index_vec) * dt; 
credit_spread_half_life_vec = 0.5/(1-prob_cutoff) * (credit_spread_recovery_index_vec - crises_starting_index_vec) * dt;
GDP_half_life_vec =  0.5/(1-prob_cutoff) *(GDP_recovery_index_vec - crises_starting_index_vec) * dt;
bank_equity_half_life_vec =  0.5/(1-prob_cutoff) * (bank_equity_recovery_index_vec - crises_starting_index_vec) * dt;

output_change_peak_to_trough = GDP_sim( crises_trough_index_vec ) ./ GDP_sim(crises_starting_index_vec) - 1;
output_change_1_years = GDP_sim(min(crises_starting_index_vec+0.5*N_interval,N_sim)) ./ GDP_sim(crises_starting_index_vec) - 1;
output_change_3_years = GDP_sim(min(crises_starting_index_vec+2.5*N_interval,N_sim)) ./ GDP_sim(crises_starting_index_vec) - 1;   % Note: the 3 year change is the average change, proxied by the mid-of-year month.
range_output_decline = quantile( output_change_1_years, 0.9 ) - quantile( output_change_1_years, 0.1 ) ;
credit_spread_half_life = mean(credit_spread_half_life_vec, 'omitnan');
bank_equity_half_life = mean(bank_equity_half_life_vec, 'omitnan');
GDP_half_life = mean(GDP_half_life_vec, 'omitnan');


%% Are credit spreads before crises too low?
% def_crisis_years = 3;
% immediately_post_crises = false * zeros( numel(1:N_interval:N_sim) , 1);
% immediately_post_distress = false * zeros( numel(1:N_interval:N_sim) , 1);
% before_crises_year_index = false * zeros( numel(1:N_interval:N_sim) , 1);
% for(i = 1:numel(crises_starting_index_yearly_vec ))
%     start_year = crises_starting_index_yearly_vec(i);
%     immediately_post_crises((start_year+1):min((start_year+def_crisis_years),length(immediately_post_crises))) = true;
%     before_crises_year_index(max(1,(start_year-def_crisis_years)):(start_year-1)) = true;
% end
% distress_years_index = find(is_distress_yearly_vec);
% for(i = 1:numel(distress_years_index ))
%     start_year = distress_years_index(i);
%     immediately_post_distress((start_year+1):min((start_year+def_crisis_years),length(immediately_post_distress))) = true;
% end
before_crises_year_index = zeros( numel(1:N_interval:N_sim) , 1);
before_crises_year_index( max(crises_starting_index_yearly_vec-1,1) ) = 1;

if(iter==1 && plot_figure)
h = figure('position', [0, 0, 600, 400]); 
subplot(1,2,1); hist(100*output_change_peak_to_trough, 7); title('Peak to Trough GDP changes');
subplot(1,2,2); hist(100*output_change_3years_in_crises_vec, 7); title('3 Year GDP Growth');
saveTightFigure(h,'Figures/Distribution of GDP In Crises')
end
% post-crisis years 
post_crisis_years = zeros( N_sim, 1 );   last_crisis_index=1;  
for( i = 1:N_sim )
    if( is_crisis_vec(i)==1 )
        last_crisis_index=i;
    end
    post_crisis_years(i) = dt * (i-last_crisis_index);
end
credit_spread_yearly = credit_spread_sim(1:N_interval:end);
credit_spread_changes_in_1_year = [ credit_spread_yearly(2:end) - credit_spread_yearly(1:(end-1)); 0  ];

crisis_prev_year = [0; is_crisis_yearly_vec(1:end-1)];
crises_post_5Y =  movsum(crisis_prev_year,  [ 4, 0 ] )>0;  % Important: the movsum function contains the summation of the current position. 
severe_crises = is_crisis_yearly_vec & (output_growth_yearly_vec<median(output_change_in_crises_vec));
mild_crises = is_crisis_yearly_vec & (output_growth_yearly_vec>median(output_change_in_crises_vec));
tbl_pre_crisis = table( before_crises_year_index, is_crisis_yearly_vec, credit_spread_yearly, severe_crises,  mild_crises,   output_growth_yearly_vec, output_growth_3year_vec,  credit_spread_changes_in_1_year,  credit_spread_changes_in_1_year.*is_crisis_yearly_vec, crises_post_5Y ,'VariableNames', ...
    { 'before_crises', 'crisis',  'credit_spread', 'severe_crises', 'mild_crises', 'output_growth', 'output_growth_in_3_years',  'credit_spread_change',  'spread_and_crisis', 'crises_post_5Y'  });
reg_spread_too_low = fitlm( tbl_pre_crisis,  'credit_spread ~ before_crises' );
reg_spread_too_low_with_control = fitlm( tbl_pre_crisis,  'credit_spread ~ before_crises + crises_post_5Y' );
reg_spread_too_low_severe = fitlm( tbl_pre_crisis,  'credit_spread ~ severe_crises + crises_post_5Y' );
reg_spread_too_low_mild = fitlm( tbl_pre_crisis,  'credit_spread ~ mild_crises + crises_post_5Y' );
reg_spread_spike =  fitlm( tbl_pre_crisis,  'output_growth_in_3_years ~ spread_and_crisis + crises_post_5Y' );



%% Conditional Regressions: Output Growth (yearly) ~ Credit Spreads
credit_spread_before_crises = credit_spread_sim( max(crises_starting_index_vec-12,1) );
credit_before_crises = bank_total_credit_sim( max(crises_starting_index_vec-12,1) );  
credit_spread_change_during_crises = credit_spread_sim( min(crises_starting_index_vec+2, N_sim) )  -  credit_spread_sim( max(crises_starting_index_vec,1) );
tbl = table( output_change_in_crises_vec, credit_spread_change_during_crises,  credit_spread_before_crises, credit_before_crises,  'VariableNames', { 'output_change',  'credit_spread_change_during_crises','credit_spread_before_crises', 'credit_before_crises' });
reg_conditional_output_and_credit = fitlm( tbl,  'output_change ~ credit_before_crises' );
reg_conditional_output_and_pre_crisis_spread = fitlm( tbl,  'output_change ~ credit_spread_before_crises' );
reg_conditional_output_and_spread_spike = fitlm( tbl,  'output_change ~ credit_spread_change_during_crises' );



%% Regressions related to equity crashes. 
pre_crash = [is_equity_crash_vec(2:end); 0];
equity_crash_next_month =  [is_equity_crash_vec(2:end); 0];
equity_crash_next_year = movsum(equity_crash_next_month,  [0,11] ) >0 ;  % Note: movsum contains the current data.
equity_crash_prev_month =  [0; is_equity_crash_vec(1:(end-1))];
post_crash_5Y = movsum(equity_crash_prev_month,  [ 12*5-1, 0 ] )>0; 
tbl = table( credit_spread_sim,  pre_crash,  post_crash_5Y, bank_credit_ratio_sim, equity_crash_next_month, equity_crash_next_year, 'VariableNames', { 'credit_spread',  'pre_crash','post_crash_5Y', 'bank_credit_GDP', 'equity_crash_next_month', 'equity_crash_next_year' });
reg_spread_pre_crash = fitlm( tbl,  'credit_spread ~ pre_crash + post_crash_5Y' );
reg_bank_credit_predicting_crash = glmfit(bank_credit_ratio_sim,  equity_crash_next_year,'binomial','Link','probit');


%% Unconditional Regressions: Output Growth (yearly) ...
% future output and current bank credit
lag=1;  
X = bank_credit_ratio_sim(1:N_interval:end); 
X1 = credit_spread_sim(1:N_interval:end); 
bank_equity_yearly = bank_equity_sim(1:N_interval:end);
X2 =  [ bank_equity_yearly(2:end)./ bank_equity_yearly(1:(end-1))-1; 0  ] ; 
X3 = current_sim.w(1:N_interval:end); 
X4 = bank_leverage_sim(1:N_interval:end);
y = double(is_crisis_yearly_vec)  ; 
y_5year = double( movsum( is_crisis_yearly_vec, [0, 4] )>0 );   % convert logical values to numeric values

output_growth_2years =  [GDP_sim_yearly(3:end)./GDP_sim_yearly(1:(end-2)) - 1; zeros(2,1)] ; 
output_growth_3years =  [GDP_sim_yearly(4:end)./GDP_sim_yearly(1:(end-3)) - 1; zeros(3,1)]; 
output_growth_4years =  [GDP_sim_yearly(5:end)./GDP_sim_yearly(1:(end-4)) - 1; zeros(4,1)]; 
output_growth_5years =  [GDP_sim_yearly(6:end)./GDP_sim_yearly(1:(end-5)) - 1; zeros(5,1)]; 

high_froth =  movmean(X1, [5,0]) < median(movmean(X1, [5,0]));  % the median of the last five years.
high_credit = X >= quantile(X,0.92);    % This follows the same definition as Krishnamurthy and Muir. (in their paper, the cutoff is 92%)
high_credit_growth = movmean(bank_credit_growth_sim_yearly, [4, 0]) > quantile(movmean(bank_credit_growth_sim_yearly, [4, 0]), 0.67);
tbl = table( X(1:(end-lag)), bank_credit_growth_sim_yearly(1:(end-lag)),  X1(1:(end-lag)),  X1(1:(end-lag))+1,  X2(1:(end-lag)), X2((1+lag):end), X3(1:(end-lag)),  X4(1:(end-lag)), high_froth(1:(end-lag)),  high_credit(1:(end-lag)), high_credit_growth(1:(end-lag)),   high_froth(1:(end-lag)).*high_credit(1:(end-lag)), y((1+lag):end), dN_yearly((1+lag):end) ,  y_5year((1+lag):end), output_growth_yearly_vec((1+lag):end),  output_growth_2years((1+lag):end),  output_growth_3years((1+lag):end), output_growth_4years((1+lag):end), output_growth_5years((1+lag):end),  is_recession_yearly_vec((1+lag):end), is_non_fin_recession_yearly_vec((1+lag):end),   ...
    'VariableNames', { 'bank_credit', 'bank_credit_growth', 'credit_spread', 'credit_spread_positive' ,'bank_equity_return', 'bank_equity_return_next_year' ,'bank_wealth_ratio',  'bank_leverage',  'high_froth', 'high_credit',  'high_credit_growth',    'high_froth_and_credit'  ,'crisis_later_1year', 'distress_later',  'crisis_later_5year'   ,'output_growth_next_year',  'output_growth_next_2years', 'output_growth_next_3years', 'output_growth_next_4years', 'output_growth_next_5years', 'recession' , 'non_fin_recession'}  );

% plot( X1, output_growth_yearly_vec , 'o',  [ credit_spread_cutoff, credit_spread_cutoff], [ min(output_growth_yearly_vec), max(output_growth_yearly_vec) ], '-.'  );
% h = lsline;  h.Color='black'; h.LineWidth=2;
% saveTightFigure( '../Figures/spread_and_future_output.pdf' )

% Predictions on output disaster. High froth and high credit. 
% reg_output_disaster_and_credit_spread= fitlm( tbl , 'crisis_later_1year ~  credit_spread')  ;
regs_output_growth_and_high_froth_and_credit.next1year = fitlm( tbl , 'output_growth_next_year ~  high_froth_and_credit');
regs_output_growth_and_high_froth_and_credit.next2year =fitlm( tbl , 'output_growth_next_2years ~  high_froth_and_credit') ;
regs_output_growth_and_high_froth_and_credit.next3year =fitlm( tbl , 'output_growth_next_3years ~  high_froth_and_credit') ;
regs_output_growth_and_high_froth_and_credit.next4year =fitlm( tbl , 'output_growth_next_4years ~  high_froth_and_credit') ;
regs_output_growth_and_high_froth_and_credit.next5year =fitlm( tbl , 'output_growth_next_5years ~  high_froth_and_credit') ;

% Froth regressions
try
    reg_output_disaster_and_froth =  fitlm(  tbl , 'crisis_later_5year ~ high_froth ');
    reg_output_disaster_and_high_credit= fitlm( tbl, 'crisis_later_5year ~  high_credit') ;
    reg_recession_and_froth = fitlm(  tbl , 'recession ~ high_froth ');
    reg_non_fin_recession_and_froth = fitlm(  tbl , 'non_fin_recession ~ high_froth ');
    reg_output_disaster_and_bank_leverage= fitlm( tbl , 'crisis_later_1year ~  bank_leverage');
    reg_output_disaster_and_bank_equity= fitlm( tbl , 'crisis_later_1year ~  bank_equity_return')  ;
catch
    reg_output_disaster_and_froth.Coefficients.Estimate = [nan;nan]; 
    reg_output_disaster_and_high_credit = 0; 
    reg_recession_and_froth = 0; 
    reg_non_fin_recession_and_froth = 0; 
    reg_output_disaster_and_bank_leverage = 0;
    reg_output_disaster_and_bank_equity = 0;
end



% reg_output_disaster_and_high_credit =  fitlm(  tbl , 'crisis_later_1year ~ high_credit ');
% reg_output_disaster_and_froth_bank_credit = fitlm( tbl, 'crisis_later_1year ~  high_froth_and_credit')  ;
% reg_output_disaster_and_froth_and_bank_credit= fitlm( tbl, 'crisis_later_1year ~  high_froth + high_credit + high_credit*high_froth ')  ;

% dependent variable: output growth
reg_output_growth_and_bank_credit= fitlm(  tbl, 'output_growth_next_year ~ bank_credit_growth')  ;
reg_output_growth_and_credit_spread =  fitlm(  tbl , 'output_growth_next_year ~ credit_spread');
reg_output_growth_and_bank_equity =  fitlm(  tbl , 'output_growth_next_year ~ bank_equity_return');
reg_output_growth_and_froth = fitlm(  tbl , 'output_growth_next_year ~ high_froth ');

% credit spread and bank credits
reg_credit_spread_and_credit = fitlm(tbl, 'credit_spread  ~ bank_credit ' )  ;




%% ****************** Bank Equity Returns and Bank Credit  ****************** 

% Bank credit ratio and future bank equity returns. 
bank_credit_ratio_yearly= bank_credit_ratio_sim( 1:N_interval:end);
X = normalize(bank_credit_ratio_yearly);
bank_equity_growth = [bank_equity_yearly((lag+1):end) ./ bank_equity_yearly(1:(end-lag)) - 1;  0] ;
bank_equity_quarterly_growth = [ bank_equity_sim( 4:3:end )  ./ bank_equity_sim( 1:3:(end-3) ) - 1; 0 ];
bank_credit_growth_last_year = [0; bank_credit_growth_sim_yearly(1:(end-1))];
y = bank_equity_growth;
y2 = [  GDP_sim_yearly((lag+1):end) ./ GDP_sim_yearly(1:(end-lag)) - 1; zeros(lag,1)] ;    % GDP growth

% IMPORTANT: Here we use -30% as the cutoff.
def_crash =-0.3; % Note that the definition of Xiong and Baron is quarterly equity decline of more than -30% 
bank_equity_crash_in_1_year = [bank_equity_quarterly_growth(5:4:end-3) < def_crash  |  bank_equity_quarterly_growth(6:4:end-2) < def_crash |  bank_equity_quarterly_growth(7:4:end-1) < def_crash |  bank_equity_quarterly_growth(8:4:end) < def_crash ;   0]  ;   % this is the same definition as Xiong and Baron

% def_crash =-0.7;   % Note: decline in a single year. This is the total decline if the quarterly decline is -30%. 
severe_bank_equity_crash_in_1_year = bank_equity_growth < -0.7; 

bank_equity_crash_in_2_year = [  bank_equity_crash_in_1_year(2:end); 0 ]  | bank_equity_crash_in_1_year;   % an equity crash in the next 2 years
bank_equity_crash_in_3_year = [  bank_equity_crash_in_1_year(3:end); 0; 0 ] | bank_equity_crash_in_2_year; 

index_large_expansion = bank_credit_ratio_yearly >=  quantile(bank_credit_ratio_yearly, 0.95);

bank_equity_return_in_1_year = [ equity_yearly((lag+1):end) ./ equity_yearly(1:(end-lag)) - 1;   zeros(lag,1) ] + par.eta + par.rho;   % equity growth
bank_equity_return_in_2_year =  [bank_equity_yearly(3:end) ./ bank_equity_yearly(1:(end-2)) - 1;  zeros(2,1) ] + par.eta + par.rho;
bank_equity_return_in_3_year =  [bank_equity_yearly(4:end) ./ bank_equity_yearly(1:(end-3)) - 1;  zeros(3,1) ] + par.eta + par.rho;


% scatter( bank_credit_ratio_yearly,  bank_equity_return_in_1_year );  xlabel('bank credit ratio'); ylabel('bank equity return in the next year');
% scatter( bank_credit_growth_last_year,  bank_equity_return_in_1_year );  xlabel('bank credit growth in the last year'); ylabel('bank equity return in the next year');

% Note: bank credit growth is not working in our setting. 
% Note: immediately_post_distress include years immediately after the last
% crisis. This is to isolate the pre-crisis behavior of the economy.
tbl = table( X, bank_equity_return_in_1_year,  severe_bank_equity_crash_in_1_year, bank_equity_crash_in_1_year, bank_equity_crash_in_2_year, bank_equity_crash_in_3_year, bank_equity_return_in_2_year, bank_equity_return_in_3_year, index_large_expansion,  'VariableNames', { 'normalized_bank_credit_ratio_yearly',  'bank_equity_return_next_year','severe_bank_equity_crash_in_1_year',  'bank_equity_crash_in_1_year', 'bank_equity_crash_in_2_year', 'bank_equity_crash_in_3_year', 'bank_equity_return_in_2_year',  'bank_equity_return_in_3_year', 'index_large_expansion' });
% Part 1: Regress bank equity crash probability on normalized bank credit ratio.  Match Table III of Xiong and Baron 2017.
reg_severe_bank_equity_crash_in_1_year_and_credit =  fitglm( tbl,  'severe_bank_equity_crash_in_1_year ~ normalized_bank_credit_ratio_yearly  ' , 'Distribution','binomial');

reg_bank_equity_crash_1Y_and_credit =  fitglm( tbl,  'bank_equity_crash_in_1_year ~ normalized_bank_credit_ratio_yearly  ' , 'Distribution','binomial');
reg_bank_equity_crash_2Y_and_credit =  fitglm( tbl,  'bank_equity_crash_in_2_year ~ normalized_bank_credit_ratio_yearly  ' , 'Distribution','binomial');
reg_bank_equity_crash_3Y_and_credit =  fitglm( tbl,  'bank_equity_crash_in_3_year ~ normalized_bank_credit_ratio_yearly  ' , 'Distribution','binomial');
% Part 2: Regress bank equity return on bank credit expansions.  Match Table III of Xiong and Baron 2017.
reg_bank_equity_and_credit = fitlm( tbl,  'bank_equity_return_next_year ~ normalized_bank_credit_ratio_yearly  ' );
reg_bank_equity_2Y_and_large_credit = fitlm( tbl,  'bank_equity_return_in_2_year ~ normalized_bank_credit_ratio_yearly  ' );
reg_bank_equity_3Y_and_large_credit = fitlm( tbl,  'bank_equity_return_in_3_year ~ normalized_bank_credit_ratio_yearly  ' );


% Equity sensitivity
sigma_X = std(bank_credit_ratio_yearly);
equity_sensitivity_to_credit = ([ X , ones(size(X)) ] \ y)  *  sigma_X; 
% equity_sensitivity_to_credit_with_crises = find_coef( X(sub_index), y(sub_index) );  equity_sensitivity_to_credit_with_crises = equity_sensitivity_to_credit_with_crises(1);
% equity_sensitivity_to_credit_without_crises = find_coef( X(~sub_index), y(~sub_index) );  equity_sensitivity_to_credit_without_crises = equity_sensitivity_to_credit_without_crises(1);
equity_sensitivity_to_credit = equity_sensitivity_to_credit(1);

sub_index = before_crises_year_index(1:(end-lag))==1;
if(iter==1&& plot_figure)
% Plot the bank equity response
h = figure('position', [0, 0, 600, 400]); 
subplot(1,3,1);plot(  X(sub_index & y<0) ,  y(sub_index & y<0), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red';  h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('With Crises'); 
subplot(1,3,2);plot(  X(~sub_index) ,  y(~sub_index), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('Without Crises');
subplot(1,3,3);plot(  X ,  y, 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('Aggregate');
saveTightFigure(h,'Figures/bank_equity_returns_and_bank_credit.pdf')

% Plot the output response. 
h = figure('position', [0, 0, 600, 400]); 
subplot(1,3,1);plot(  X(sub_index) ,  y2(sub_index), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red';  h2.LineStyle='-.'; ylabel('output growth in 3 years'); title('With Crises'); 
subplot(1,3,2);plot(  X(~sub_index) ,  y2(~sub_index), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('output growth in 3 years'); title('Without Crises');
subplot(1,3,3);plot(  X ,  y2, 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('output growth in 3 years'); title('Aggregate');
saveTightFigure(h,'Figures/output_growth_and_bank_credit.pdf')
end

% Replicate Xiong and Baron Table III. The equity risk premium turns
% negative in their model when the credit expansions are large. 
% mean( bank_equity_return_sim_yearly(bank_credit_ratio_yearly>5) )
% mean( bank_equity_return_sim(bank_credit_ratio_sim<=5 & bank_credit_ratio_sim>3) )*12
% mean( bank_equity_return_sim(bank_credit_ratio_sim<=3 & bank_credit_ratio_sim>1) )*12
% mean( bank_equity_return_sim(bank_credit_ratio_sim<=1) )*12
avg_bank_equity_return_conditional_on_medium_credit = mean( bank_equity_return_sim(bank_credit_ratio_sim>=quantile(bank_credit_ratio_sim,0.2) & bank_credit_ratio_sim<quantile(bank_credit_ratio_sim,0.9) ) + par.eta  ) ;
avg_bank_equity_return_conditional_on_low_credit = mean( bank_equity_return_sim(bank_credit_ratio_sim<quantile(bank_credit_ratio_sim,0.2) ) + par.eta  ) ;



%% Bank Equity Returns and Bank Leverage
% What about bank leverage?
X= bank_leverage_sim( 1:N_interval:end );  X = X/std(X);
tbl = table( X, y, 'VariableNames', { 'bank_leverage_standardized',  'bank_equity_return_in_2_years' });
reg_bank_equity_and_bank_leverage = fitlm( tbl,  'bank_equity_return_in_2_years ~ bank_leverage_standardized  ' );
bank_equity_to_bank_leverage_sensitivity = reg_bank_equity_and_bank_leverage.Coefficients(2,1).Estimate;

sub_index = before_crises_year_index(1:(end-lag))==1;
if(iter==1&& plot_figure)
h = figure('position', [0, 0, 600, 400]); 
subplot(1,3,1);plot(  X(sub_index) ,  y(sub_index), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red';  h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('With Crises'); 
subplot(1,3,2);plot(  X(~sub_index) ,  y(~sub_index), 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('Without Crises');
subplot(1,3,3);plot(  X ,  y, 'o', 'LineWidth', 1 ); h2 = lsline; h2.LineWidth = 2;  h2.Color='red'; h2.LineStyle='-.'; ylabel('bank equity returns in 3 years'); title('Aggregate');
saveTightFigure(h,'Figures/bank_equity_returns_and_bank_leverage.pdf')
end

%% ----------- Key Features of the Model ----------- 
% Feature 1: Calm Before the Storm
index_pre_crises =  [ max(1,crises_starting_index_vec-1) ];
index_after_crises = [];
for(i = 12:(20*12))
    index_after_crises = [ index_after_crises;  min(crises_starting_index_vec+i, N_sim) ];  % exclusion of the 
end
index_after_crises = unique(index_after_crises);

vol_sim = current_sim.vol;    vol_premia_sim = current_sim.vol_premia;
vol_before_crises = vol_sim(index_pre_crises);
vol_premia_before_crises = vol_premia_sim(index_pre_crises);
vol_ratio = mean(vol_before_crises,'omitnan') / mean(vol_sim(index_after_crises),'omitnan');
vol_premia_ratio = mean(vol_premia_before_crises,'omitnan') / mean(vol_premia_sim(index_after_crises),'omitnan');
feature1.vol_ratio = vol_ratio;
feature1.vol_premia_ratio = vol_premia_ratio;

% Feature 2: Predictability of Severe Crises
lag_month = 12;
credit_growth = [ mean(bank_credit_ratio_sim(2:end)./bank_credit_ratio_sim(1:(end-1))-1);  bank_credit_ratio_sim(2:end)./bank_credit_ratio_sim(1:(end-1))-1 ];
output_growth = [ mean(current_sim.GDP(2:end)./current_sim.GDP(1:(end-1))-1);  current_sim.GDP(2:end)./current_sim.GDP(1:(end-1))-1 ];
y = double(is_crisis_vec)  ;   % Note: need to convert the logical values into numerical values.
y1 = double(dN_effective==1) ;
tbl = table( bank_credit_ratio_sim(1:(end-lag_month)), credit_growth(1:(end-lag_month)), credit_spread_sim(1:(end-lag_month)),  credit_spread_sim((1+lag_month):end) , bank_equity_return_sim(1:(end-lag_month)), current_sim.w(1:(end-lag_month)), current_sim.xK(1:(end-lag_month)), y((1+lag_month):end), y1((1+lag_month):end), 'VariableNames', { 'bank_credit', 'bank_credit_growth', 'credit_spread', 'credit_spread_future', 'bank_equity_return', 'bank_wealth_ratio', 'bank_leverage', 'crisis_later_1year', 'crisis_indicator'});

% [coef, dev, stats ] = mnrfit(   tbl.bank_leverage,  ordinal(tbl.crisis_indicator, {'1', '0'}, [true, false] ) , 'model', 'ordinal'   )
reg_severe_crises_and_bank_leverage= fitlm( tbl , 'crisis_later_1year ~  bank_leverage')  ;
reg_crises_and_bank_leverage= fitlm( tbl , 'crisis_indicator ~  bank_leverage')  ;
reg_bank_credit_growth_and_bank_leverage = fitlm( tbl , 'bank_credit_growth ~  bank_leverage')  ;
reg_bank_credit_spread_future_and_bank_leverage = fitlm( tbl , 'credit_spread_future ~  bank_leverage') ;

feature2.reg_severe_crises_and_bank_leverage=reg_severe_crises_and_bank_leverage;
feature2.reg_crises_and_bank_leverage=reg_crises_and_bank_leverage;
feature2.reg_bank_credit_growth_and_bank_leverage=reg_bank_credit_growth_and_bank_leverage;
feature2.reg_bank_credit_spread_future_and_bank_leverage=reg_bank_credit_spread_future_and_bank_leverage;

% Feature 3: Pre-crisis credit spread too-low
tbl = table( current_sim.GDP( min(crises_starting_index_vec+1, N_sim) )./current_sim.GDP(crises_starting_index_vec)-1, ...
    credit_spread_sim(min(crises_starting_index_vec+1, N_sim)) - mean(credit_spread_sim) , ...
    credit_spread_sim(index_pre_crises), current_sim.xK(index_pre_crises),  ....
    'VariableNames', { 'output_change',  'credit_spread_change_during_crises','credit_spread_before_crises', 'bank_equity_before_crises' });
reg_pre_crisis_low_spread = fitlm( tbl , 'output_change ~  credit_spread_change_during_crises + credit_spread_before_crises');
feature3.reg_pre_crisis_low_spread = reg_pre_crisis_low_spread;

% Feature 4: bank equity returns are too low. 
bank_equity_ratio = mean(bank_equity_return_sim(index_pre_crises),'omitnan') / mean(bank_equity_return_sim(index_after_crises),'omitnan') - 1;
feature4.bank_equity_ratio = bank_equity_ratio;

% Feature 5: bond returns predictability
w_dim = kron(ones(par.N_lambda,1), par.w_grid );
lambda_dim = kron(par.lambda_grid, ones(par.N_w, 1) );
fit_a_function = @(values_matrix)  scatteredInterpolant(w_dim,    lambda_dim ,  reshape( values_matrix , par.N_w*par.N_lambda, 1 ),   'linear',   'linear')  ;   % note: linear extrapolation at the boundary
excess_bond_return_matrix = model_sol.lambda_matrix .* 0.5 ./(1 - model_sol.kappab_survival_matrix  );
bond_return_fun = fit_a_function( excess_bond_return_matrix );
bond_excess_return_sim = bond_return_fun( current_sim.w, current_sim.lambda );

index_future = 13:N_sim;  index_current = 1:(N_sim-12);
tbl = table(  bond_excess_return_sim(index_future), bank_credit_ratio_sim(index_current), bank_leverage_sim(index_current), bank_equity_return_sim(index_current)  , 'VariableNames', {'bond_excess_return', 'bank_credit', 'bank_leverage', 'bank_equity_return' });
reg_bond_return_prediction = fitlm( tbl , 'bond_excess_return ~  bank_credit + bank_leverage + bank_equity_return ');
feature5.reg_bond_return_prediction = reg_bond_return_prediction;



%% Store Key Results
scenario = current_sim;
scenario.w_avg = mean(current_sim.w);

% Crises index 
scenario.crises_starting_index_vec = crises_starting_index_vec;
scenario.crises_starting_index_yearly_vec = crises_starting_index_yearly_vec;
scenario.is_crisis_yearly_vec = is_crisis_yearly_vec;
scenario.is_distress_yearly_vec = is_distress_yearly_vec;
scenario.is_crisis_vec = is_crisis_vec;

% Frequency of crisese before eliminating clustered ones. 
scenario.frequency_non_sparse_crisis = frequency_non_sparse_crisis;
scenario.frequency_non_sparse_equity_crashes = frequency_non_sparse_equity_crashes;


% Store key moments. 
scenario.frequency_crises_simulated = frequency_crises_simulated;
scenario.frequency_distress_underlying = frequency_distress_underlying;
scenario.frequency_distress_simulated = frequency_distress_simulated;

scenario.bank_equity_change_in_crises_avg = bank_equity_change_in_crises_avg;
scenario.bank_equity_yearly_return_avg = bank_equity_yearly_return_avg;
scenario.bank_equity_yearly_return_vol = bank_equity_yearly_return_vol;
scenario.bank_equity_excess_return_sim_yearly = bank_equity_excess_return_sim_yearly;

scenario.output_change_in_crises_avg = output_change_in_crises_avg;
scenario.range_output_decline = range_output_decline;
scenario.output_growth_avg = output_growth_avg;
scenario.output_growth_vol = output_growth_vol;
scenario.output_growth_monthly_vec = output_growth_monthly_vec;
scenario.bank_credit_growth_monthly = bank_credit_growth_monthly;
scenario.is_equity_crash_vec = is_equity_crash_vec;

scenario.bank_credit_GDP_ratio_sim = bank_credit_ratio_sim;  % note: after the change on Nov 15, bank_creidt_ratio is measuring the credit/GDP ratio.  

scenario.credit_spread_change_in_crises_avg = credit_spread_change_in_crises_avg;

% capital excess return, to proxy for the risk premium.
p = current_sim.p ;    K=current_sim.K;    pK = p.*K;
% scenario.capital_excess_return = [ current_sim.p(13:end)./current_sim.p(1:(end-12))- 1; zeros(12,1)  ]; % - current_sim.investment_rate./current_sim.p + par.AH./current_sim.p  - current_sim.rf;
flow_vec = (- current_sim.investment_rate./current_sim.p + par.AH./current_sim.p  - current_sim.rf) * dt;
flow_vec = movsum(  flow_vec, [0,12]  );  % Note: move sum the next 12 elements, not BEFORE.  This detail is very important. 
%****** Important: capital excess return from the model.   Need to take care of the alpha component. 
% capital_excess_return = [ pK(13:end)./pK(1:(end-12)) - 1 - par.alpha * movsum(dN_effective .* (current_sim.xK-1) ./ (current_sim.xK) , [0,12], 'Endpoints','discard' )  ; zeros(12,1)  ] + flow_vec; %- current_sim.investment_rate./current_sim.p + par.AH./current_sim.p  - current_sim.rf;
capital_excess_return = [ pK(13:end)./pK(1:(end-12)) - 1; zeros(12,1)  ] + flow_vec; 
capital_excess_return = capital_excess_return - mean(capital_excess_return) + mean(current_sim.capital_risk_premium_vec)*12;
scenario.capital_excess_return = capital_excess_return;
scenario.capital_excess_return_yearly = scenario.capital_excess_return(1:12:end);

high_credit_quantile=0.8;

tbl = table(  bank_credit_ratio_yearly(1:(end-1)), scenario.capital_excess_return_yearly(2:end), scenario.bank_equity_excess_return_sim_yearly(2:end), 'VariableNames', { 'bank_credit', 'capital_excess_return','bank_equity_exess_return'});
reg_capital_excess_return_and_bank_credit= fitlm( tbl , 'capital_excess_return ~  bank_credit')  ; 
scenario.reg_capital_excess_return_and_bank_credit = reg_capital_excess_return_and_bank_credit;
scenario.capital_excess_return_on_high_credit =  mean(scenario.capital_excess_return_yearly( bank_credit_ratio_yearly >  quantile( bank_credit_ratio_yearly, high_credit_quantile )  )) ;

reg_bank_equity_excess_return_and_bank_credit= fitlm( tbl , 'bank_equity_exess_return ~  bank_credit')  ; 
scenario.reg_bank_equity_excess_return_and_bank_credit = reg_bank_equity_excess_return_and_bank_credit;
scenario.bank_equity_excess_return_on_high_credit =  mean(scenario.bank_equity_excess_return_sim_yearly( bank_credit_ratio_yearly >  quantile( bank_credit_ratio_yearly, high_credit_quantile )  )) ;


% volatility of output growth during a crisis shock. 
scenario.output_growth_volatility_in_shocks = std(  output_change_1_years   );

% recovery speeds
scenario.credit_spread_half_life = credit_spread_half_life;
scenario.GDP_half_life = GDP_half_life;
scenario.bank_equity_half_life = bank_equity_half_life;

% Yearly values
scenario.credit_spread_yearly = credit_spread_sim(1:N_interval:end); 
scenario.output_growth_yearly = output_growth_yearly_vec;   % Note: at year t, this measures the output growth from t to t+1. In otherwords, we meausre previsible output growth.
scenario.bank_credit_ratio_yearly = bank_credit_ratio_sim(1:12:end);
scenario.bank_credit_growth_yearly = bank_credit_growth_sim_yearly;

% post-crisis years 
scenario.post_crisis_years = post_crisis_years;

% regressions
scenario.reg_output_and_bank_credit = reg_output_growth_and_bank_credit;
scenario.reg_output_disaster_and_bank_credit = reg_output_disaster_and_high_credit;
scenario.reg_output_disaster_and_bank_equity = reg_output_disaster_and_bank_equity;
scenario.reg_output_disaster_and_bank_leverage = reg_output_disaster_and_bank_leverage;

scenario.reg_bank_equity_and_credit = reg_bank_equity_and_credit;
scenario.reg_credit_spread_and_credit = reg_credit_spread_and_credit;
scenario.reg_conditional_output_and_credit = reg_conditional_output_and_credit;
scenario.reg_conditional_output_and_pre_crisis_spread=reg_conditional_output_and_pre_crisis_spread;
scenario.reg_conditional_output_and_spread_spike=reg_conditional_output_and_spread_spike;
scenario.reg_spread_spike = reg_spread_spike;   % this is a regerssion on the interaction between credit spread and the indicator of crisis.
scenario.regs_output_growth_and_high_froth_and_credit = regs_output_growth_and_high_froth_and_credit;

% bank equity crahses. 
scenario.reg_spread_pre_crash = reg_spread_pre_crash;
scenario.reg_bank_credit_predicting_crash = reg_bank_credit_predicting_crash; 

% bank equity related. 
scenario.reg_severe_bank_equity_crash_in_1_year_and_credit = reg_severe_bank_equity_crash_in_1_year_and_credit;
scenario.reg_bank_equity_crash_1Y_and_credit = reg_bank_equity_crash_1Y_and_credit;
scenario.reg_bank_equity_crash_2Y_and_credit = reg_bank_equity_crash_2Y_and_credit;
scenario.reg_bank_equity_crash_3Y_and_credit = reg_bank_equity_crash_3Y_and_credit;

scenario.reg_bank_equity_1Y_and_large_credit = reg_bank_equity_and_credit;
scenario.reg_bank_equity_2Y_and_large_credit = reg_bank_equity_2Y_and_large_credit;
scenario.reg_bank_equity_3Y_and_large_credit = reg_bank_equity_3Y_and_large_credit;
scenario.bank_equity_growth = bank_equity_growth;
scenario.avg_bank_equity_return_conditional_on_medium_credit = avg_bank_equity_return_conditional_on_medium_credit;
scenario.avg_bank_equity_return_conditional_on_low_credit=avg_bank_equity_return_conditional_on_low_credit;

% spread too low? 
scenario.reg_spread_too_low = reg_spread_too_low;
scenario.reg_spread_too_low_mild = reg_spread_too_low_mild;
scenario.reg_spread_too_low_severe = reg_spread_too_low_severe;
scenario.reg_spread_too_low_with_control  = reg_spread_too_low_with_control;

% Core regressions for distinguishing rational and irrational
% scenario.reg_output_disaster_and_credit_spread = reg_output_disaster_and_credit_spread;
scenario.reg_output_growth_and_credit_spread = reg_output_growth_and_credit_spread;  % unconditional regressions of output growth on credit spreads
scenario.reg_output_disaster_and_froth = reg_output_disaster_and_froth;
scenario.reg_output_growth_and_froth = reg_output_growth_and_froth;
scenario.reg_output_growth_and_bank_equity = reg_output_growth_and_bank_equity;

% others
scenario.equity_sensitivity_to_credit = equity_sensitivity_to_credit;  % note: standardized bank credit. So if the coefficient is -1%, then it says if bank credit increases by one standard deviation, bank equity in the next 2 years go down by 1%.
scenario.bank_equity_to_bank_leverage_sensitivity = bank_equity_to_bank_leverage_sensitivity;   % bank equity sensititivy to bank leverage, in two years.  
scenario.output_to_capital_ratio_avg = mean(current_sim.productivity./current_sim.p);
scenario.bank_lvg_avg = mean( current_sim.xK );

% store vectors
scenario.output_change_in_crises_vec = output_change_in_crises_vec;
scenario.output_change_peak_to_trough = output_change_peak_to_trough;
scenario.output_change_1_years = output_change_1_years;
scenario.output_change_3_years = output_change_3_years;

scenario.recovery_length_vec = recovery_length_vec;
scenario.peak_to_trough_vec = peak_to_trough_vec;

% features 
scenario.feature1 = feature1;
scenario.feature2 = feature2;
scenario.feature3 = feature3;
scenario.feature4 = feature4;
scenario.feature5 = feature5;

% average paths of output and equity, credit spread.
scenario.average_path = average_path;

scenario.dN_vec = sim.dN_vec;
scenario.dB_vec = sim.dB_vec;

% Sharpe ratios of capital returns from bank's perspective. 
investment_growth_return = (par.AH-current_sim.investment_rate) ./ current_sim.p ;  % Note: this is a yearly value
RK_monthly = diff(log(current_sim.p) + current_sim.logK_vec_original) +  investment_growth_return(1:(end-1))/12 ;
RK_monthly = [mean(RK_monthly); RK_monthly];  % make sure the vector length is the same. 
scenario.RK_monthly = RK_monthly;

% a table of data. 
scenario.tbl_output = table(  sim.T_vec,  current_sim.w,   current_sim.lambda,   current_sim.lambda_rational,  current_sim.K,  current_sim.p,  current_sim.GDP, current_sim.productivity  , current_sim.credit_spread, current_sim.bank_equity, current_sim.xK,  current_sim.vol,  current_sim.vol_premia ,    sim.dN_vec,  is_crisis_vec,  sim.dB_vec,  current_sim.rf,  current_sim.rK, current_sim.rd,  bank_equity_excess_yearly_return_sim,  current_sim.sharpe_ratio,  current_sim.investment_rate,  scenario.capital_excess_return, current_sim.bank_credit,   RK_monthly, current_sim.kappap,  current_sim.logK_vec_original, ...
                'VariableNames', {   'T',                   'w',                     'lambda',              'lambda_rational',                     'K',                'p',                    'GDP',       'productivity'               , 'credit_spread',                    'bank_equity',                  'xK',                   'vol_of_p',                      'dB_premia',           'dN',          'crisis'  ,         'dB',                'rf',                  'rK',                       'rd',                       'equity_excess_return',                                            'sharpe_ratio',                      'investment_rate',                          'capital_excess_return',               'bank_credit',  'rK_monthly_from_sim',  'kappap' , 'logK_vec_original'}   );

if(iter==1)
    rational_moments = scenario;
else
    behavioral_moments = scenario;
end


end




